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An analytical study of the return time distribution of extreme events for stochastic processes 
with power-law correlation has been carried on. The calculation is based on an e-expansion in the 
correlation exponent: C(t) = jrJ~ 1+E . The fixed point of the theory is associated with stretched 
exponential scaling of the distribution; analytical expressions, valid in the pre-asymptotic regime, 
have been provided. Also the permanence time distribution appears to be characterized by stretched 
exponential scaling. The conditions for application of the theory to non-Gaussian processes have 
been analyzed and the relations with the issue of return times in the case of multifractal measures 
have been discussed. 



PACS numbers: 02.50.Ey,05.40.Ca,05.45.Tp,89.75.-k 

I. INTRODUCTION 

Calculating the return time statistics of rare events 
in stochastic processes, is one of the classical problems 
in probability theory. The applications are widespread, 
one of the most direct being the determination of safety 
margins against catastrophic events like floods and earth- 
quakes. In other fields, like statistical mechanics and in- 
formation theory, the statistics of return times plays an 
important role, being intimately connected with the way 
a system looses memory of its initial conditions. 

Starting from the work of D6eblinll| , and that of Bell- 
man and Harris on Markov chains [2], one of the known 
results is that, if the system correlations decay sufficiently 
fast, the distribution of the return times of asymptoti- 
cally rare events will tend to be exponential (see [3| for 
recent references). 

More recently, motivated by the observation that a 
wide variety of experimental records present long-time 
correlations (see |4| and references therein) there has 
been growing interest in the case the stochastic process 
is power-law correlated, with exponent small enough for 
the correlation time to be infinite. An interesting result is 
that, in this case, the return time distribution of extreme 
events appears to be well fitted by a stretched exponen- 
tial, with exponent equal to the power in the correlation 
decay, rather than by a simple exponential [1, H, H, 0] • 

Actually, seen in the light of the classical work by 
Newell and Rosenblatt @ on the probability of no zero 
crossing for power-law correlated Gaussian processes, this 
result is not really surprising. If the rare event is associ- 
ated with the crossing of a high threshold by a stochas- 
tic variable, the return time probability will be the no- 
crossing probability, for an initial condition in which the 
variable is right below threshold. 

The central idea in Q is that threshold crossings will 
be less likely for processes with longer correlations. Com- 
paring with processes with known return time distribu- 
tion (typically, a superposition of an Ornstein-Uhlembeck 
process and a stochastic variable) Newell and Roseblatt 
were able to prove that the no zero-crossing probability 
for y is bounded from above and from below by stretched 



exponentials. 

In principle, the approach in Q could be extended to 
the case of a threshold different from zero, allowing to 
conclude that stretched exponential is the correct asymp- 
totic scaling of the return probability for large values of 
its argument. This leaves open, however, some important 
questions. It would be of obvious interest to tighten the 
inequalities in Q, fixing the values of the prefactors in 
the stretched exponential scaling. It would also be inter- 
esting to have some idea of how and when (and perhaps 
why) the asymptotic regime is reached, and if Gaussian- 
ity is really an essential hypothesis. 

Purpose of this paper is to present an analytical treat- 
ment of these isssues, based on renormalized perturba- 
tion theory and an e-expansion in the correlation expo- 
nent: (y(t)y(0)) ~ |i| _1+e ; this is basically an expansion 
around the transition to infinite correlation time. The 
expansion will turn out to work well also for rather large 
values of e, providing in the Gaussian case, valid approx- 
imate expressions for the return time distributions. 

Among the other things, the analysis will point out 
the dominance of transient behaviors, in any range of 
practical interest for the return times. It will also point 
out that, perhaps contrary to intuition, the return time 
distribution is less sensitive to the extreme statistics of 
the process than to its correlation structure, in particular 
that of the correlations between scales. This will allow 
extension of the results to rather generic non-Gaussian 
processes. 

This paper is organized as follows. In Section II, the 
main definitions and results are recalled. In Section III, it 
is shown how the long memory of the process is associated 
with secular behaviors in the evolution equations for the 
exit probabilities. Section IV contains the main results 
of the paper, and analytic expressions for the return time 
distribution are provided in the form of a renormalized e- 
expansion. Sections V and VI focus again on the relation 
between the statistics of the stochastic process and that 
of the return times, extending the results to the case of 
non-Gaussian processes. Section VII is devoted to the 
statistics of permanence above threshold. Section VIII 
contains the conclusions. 
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II. RETURN, PERMANENCE AND EXIT 

Let us consider a unit variance and zero mean, sta- 
tionary Gaussian process y(t), with correlation C{t) = 
(y(t)y(0)) decaying like a power-law at sufficiently long 
time separations: 

C{t) = A\t\- 1+ \ |t| > t-o, 0<e<l. (1) 

This corresponds to the energy spectrum: 

C u = J dte iuJt C(t) ~ Lu~\ (2) 

which is a common occurrence in several physical systems 
Q. Notice also that the scaling in Eq. (p} is that of 
the velocity of a superdiffusive particle: y — r, (\r(t) — 
r(0)| 2 ) ~ |£| 2 ~ e - For e < 0, in turn, the correlation time 
would be 0(tq), and (\r(t) — r(0)| 2 ) ~ \t\, like in the case 
of a Brownian particle. For e > 1, Eq. ^ would give 
the spectrum of a signal with local anomalous diffusive 
behavior (\y(t) - y(0)\ 2 ) ~ \t\~ 1+e . 

The Gaussian process y(t) could be generated numer- 
ically, e.g. by the algorithm described in (Toj |. approxi- 
mating y(t) by a superposition of independent Ornstein- 
Uhlcnbeck processes x n (t), n — 0,1, ...N. A correlation 
like the one in Eq. ([1]) could be obtained setting for the 
correlation time and variance of x n : 

t„ = t 2" and < = L__^__. ( 3 ) 

The lower line in Fig. Q]is the correlation profile obtained 
with this technique, for e = 0.5 and TV = 24. 

We identify an extreme event by the condition y > q, 
where q is a sufficiently large threshold for the process 
and introduce two occurrence times: the permanence 
time S q of the variable y above threshold, and its coun- 
terpart below threshold, the return time R q to the event, 
which coincides with the first exit time from y < q, for 
initial condition y(0) = q. The averages R q and S g , can 
be related to the event probability P(y > q) by means of 
the general relation (Kac theorem [HI]): 

S q /R q ~P(y>q). (4) 

We are interested in the occurrence time distributions 
P(S q > t) and P(R q > t). We focus first on the return 
time distribution P(R q > t) and we can write: 

P(R q >t) = P(y(r)<q,r€[0,t]\y(0) = q), (5) 

which coincides with the no-exit probability for initial 
condition y(0) = q. This probability is the integral from 
t = to t = t of the probability current across q. If y(t) 
obeyed a stochastic differential equation, the whole prob- 
lem would reduce to solution of a Fokker-Planck equation 
with absorbing boundaries at q 

In general, calculation of the current requires knowl- 
edge of the profile near q of the conditional PDF (proba- 
bility density function) p(y(t)\y(r) < q,r S [0,i];y(0) = 
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FIG. 1: Scaling of the correlation function C(i) (lower line) 
and the probability difference AP(t) = P(y(t) > q\y(0) > 
q) - P(y > q) (upper line), for e = 0.5 and N = 24. The thin 
line is t~ ' 5 . The correlation function is obtained from Eq. 
using C(t) — J2n( Xn ( t ) Xn W}'^ ^P(t) is obtained using 
this expression for C(t) and Eqs. (|1H12|I . 



q). Exponential scaling of P(R q > r) would be associated 
with a current —jP(R q > t), with 7 a constant, deter- 
mined by the limit form of p(y(t)\y{r) < q,r G [— oo,i]), 
when the initial condition time is sent to — 00. In the 
case of a long correlated stochastic process, it will ap- 
pear that this limit form is associated with zero current 
and that the approach to the limit is what generates the 
anomalous scaling of the return time distribution. 

III. THE EFFECT OF MEMORY 

The problem becomes more tractable if sampling is 
carried on at discrete times tk = kA; indicate then y^ = 
y(t k ) and define A = AA~ 1+e , so that C{t k ) = A\k\- l+t . 
Notice that for A > tq, the statistics would become itself 
A dependent, and, for large A: S q (A) ~ A; from Kac 
theorem we would then have: R q (A) ~ A/P(y > q). 
Normalizing times with respect to A, the discrete version 
of Eq. (O will read: 

P(R q >n) = P(y k <q,k = l, ...n\y > q), 

which can be expressed in terms of the conditional prob- 
abilities: 

n 

P{y k <q,k = l, ...n\y > q) = - P k (k)], (6) 

fe=i 

with P k {k) = P(y k > q\yi < q,l = I, ..k - l;y > q) the 
return probability at time t k conditioned to being a first 
return. More in general, introduce the return probability 
conditioned to no exits before a time t\ < t n : 

Pi(n) = P(y n > q\yj <q,j = 1, ...I -l;y > q), (7) 
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and the return probability conditioned to exits at times 
tfci < tk 2 < -4 P ! an d no exits before a time t\ < t^: 



P;(n|k) = Pi(n\y kt > q,i = l,...p). 



(8) 



From here, an exact recursion relation can be derived, 
describing the evolution of the return probability from 
the initial condition at n = 0. Equations (05]) provide 
us, in fact, with the following identities: 

P n +i(m) = P n {m\y n < q) 

= [1 - Pn(n)}^ 1 P n (y m > q,y n < q) 

= [1 - P^n)]- 1 ^™) - P n (y n . m > ?)] 

= [1 - P n (n)]- 1 [P n (m) - P n (n)P n (m|n)], 

which can be rearranged to give: 

P n +i(m) = Pn(m) + P n (n)[P n (m) - P n (m\n)}, 

where P n {n) = P n (n)[l — P n (n)] -1 . In a certain sense, 
this is the analog, for a discrete non-Markovian process, 
of the Fokker-Planck equation with absorbing boundary 
conditions discussed previously, and can be iterated to 
give: 

n-1 

P„(m) = P x (m) + J2 P(0[PM - PiHl)]. (9) 



The physical meaning of this equation is to quantify how 
many of the phase points above threshold at times / = 
1, ...n— 1 should be subtracted from the probability mass, 
that, without taking into account the condition of first 
return, would be above threshold at time to. 

In order to solve Eq. © , we need to know the function 
Pl(m\l), which is essentially the second return probabil- 
ity. Repeating the same steps leading to Eq. ([9|) with 
the probability P;(m|k), we see that Eq. © is the first 
item in an unclosed hierarchy of equations, whose generic 
element reads: 



P„(m|k) = P(m|k) + PG|k)[PHk) - P/(m|k0], 

1=1 

(10) 

where fl(/|k) = Pi(l\k)[l - fl(l|k)] -1 - 

A natural strategy could be, at this point, perturbation 
theory around the "ground states" P^ 0) (m) = Pi(m) and 
P^ 0) (m|k) = Pi(m|k). Notice that substituting P k (k) -> 
Pi (k) in Eq. ((6|) and sending n — > oo would lead to 
exponential scaling of the return time distribution: 



P{R q > n) ~ [1 - Pi(oo)]" ~ exp(-n/P g ), 

where R q ~ 1/P(y > <z) [compare with Eq. Q] and we 
have exploited Pi(oo) = P(y > q) > 0. 

The lowest order expressions P^ y) (m) = Pi(m) and 
P^ 0) (m|k) = Pi (to | k) can be calculated explicitly from 
PDF's in the form: 



Px{n) 



1 



(2 7 r) 1 /2 (Tl ( n ) 



exp i 



(y n - ^i(n)) 2 



(11) 



and similar expression for pi(n\k). For large g, in fact, 
the conditions y , y ki > q in Pi(k) and Pi(n|k) can be re- 
placed by ?/o, yki = q and this guarantees Gaussian statis- 
tics. The conditional mean and variance in Eq. can 
be written in the following form (see e.g. [13| . Appendix 
C): 



Mi(n|k) = gEy C-iAj, 
ff 1 (n|k)=E«C t AiC;. 



(12) 



where Cj = C(n — ki), Dij is the inverse of the matrix 
C (ki — kj) and we have defined fco = so that now i,j = 
0,1, ...p. 

Unfortunately, we are going to see that the first order 
correction ?W(m) = Yd=i A (0 [ p l M ~ Pl HO] di- 
verges for n — > oo, and the same occurs with P^ l5 (m|k). 
As expected, exponential scaling does not appear to be 
an appropriate guess for the asymptotic behavior of the 
return time distribution. 

Let us prove this. Rearrange indices so that 
Pi(0[Pi(m)-Pi(m|0] -> Pi{m-l)[P x {m)-Px{m\m-l)] 
and send n, to — * oo with to — n finite, so that memory 
of the initial condition at t = is lost. We have in this 
limit: 

oo 

P£ ) (m)= E PiMlPiM-P^]. 

l=m— n+l 

The PDF jOi(oo) associated with Pi(oo) is just the equi- 
librium PDF for y, corresponding to setting in Eq. (fTTj) 
(j, = and a 2 = 1. From Eqs. |(T2J) and fT) wc find: 



= qAl 



and o{{l) = \-A z l 



2;-2+2e 



(13) 



and, for large I, we can Taylor expand Eq. pip : 

pi(0 = (ah(0w + \{yl - i)K(0 - i))pi(oo). 

Substituting into Pi(Z) = J^ 00 pi(l)dyi and calculating the 
integral by the Laplace method, we find: 

PiM-PW = -( g 2 /2)Pi(oo)^- 1+e +0(Z- 2 - 2e ) (14) 

where Pi(oo) ~ (27r) — 1 / 2 g _1 exp(— q 2 /2); this scaling can 
be compared with the top curve in Fig. [1] Substituting 
into 53 Pi (oo) [Pi (oo) — Pi(l)], we find that the sum is 
divergent for e > 0. Notice that the leading contribution 
to the scaling in Eq. (fl"4]) . the one that causes divergence 
of P^ ) {m), is the slow decay of the conditional mean 
Mi(0 in Eq. CE3D. 

Let us prove that also P^ i:) (TO|k) diverges as the initial 
condition is sent to t — — oo. The summand in P^ 1 ' (m|k) 
[see Eq. (fTD"])] reads: 

Pi(Z|k)[Pi(m|k)-Pi(rn|kO] 

= Pi(m|k)(Pi(Z|k)/Pi(Z|k))[Pi(Z|k) - Pi(Z|mk)]. 

From Eq. (JT2J), for fe x - I > to- fa: ^i (%), /ii(Z|mk) ~ 
gC(fa — OEii Aj ~ (fa — _1+e an( i we fi n( l again 
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f\(Z|k),i\(J|mk) ~ [1 + 0(|*i - /|- 1+£ )]Pi(oo). This 
leads to the expression: 

n — 1 n 

J2 Aa|k)[Pi(m|k) - Px(m\kl)] ~ £ r 1+£ , (15) 



/=i 



which grows like n £ ; again, divergence is produced by the 
slow decay of the conditional mean 



IV. e-EXPANSION 

We have seen that a "bare" perturbation expansion of 
Eqs. (Pl fTUl) . around the zeroth order: P^ 0) (to) = Pi(m) 
and P^ 0> (m|k) = Pi(m|k) leads to infinities. Some kind 
of renormalization is necessary, however, in order to have 
a workable theory, we still need that, to lowest order, the 
equations in the hierarchy pO[) remain decoupled. This 
basically fixes the renormalization procedure. 



We separate out of the probabilities in Eqs. (|9l fT0|) a 
renormalized part P R and a remnant P N : P = P R + P N , 

and expand P R = Ez=o A ' pR ' where A o pR = p i- Thc 
first order renormalization to Eq. (J9j) reads: 



n-l 



AiP«(m) = > [AiPf (0 + P 1 (l)][P 1 (m) - P(m|0] s 



i=i 



where [Pi (to) — Pi(to|^)] scc contains the contribution 
leading to divergence of P^ 1 ' (to), and, seeking an analogy 
with quantum field theory (see also [lj]), AiP, R (Z) could 
loosely be seen as a counterterm. Analogous expressions 
are obtained for AiP R (TO,|k) and Eq. (JTUJ) . Combining 
with Pi(m), we get the lowest order equation for the 
renormalized probability: 



n-l 



P r (to) = Pi(m) + > P ; R (0[PiM - Pi(r 



i=i 



(16) 



The p-th order renormalization, can be expressed in the 
form: 

A p P^(to) = E?=i{\Pm[Pi(m) - PMOUSc 1 ' 

+P r ' p (0A p _![P(to) - P(to|0] scc }, (17) 

and we use here a superscript to indicate the order 
at which each expression is considered: P np = P x + 
Ym=i A;P R and [. . = Ef=o M- Equation 

(|17[) can be rewritten in the more concise form, which 
generalizes Eq. (fT6|) : 



P R ^(TO) = P(m) 



n-l 

E 

i=i 



(OIPM-PM 



i(p~i> 

J sec 



Again, analogous expressions hold for A p P R (TO,|k) and 
Eq. (fini). From inspection of Eqs. (TO]) and (TTTJ), we 
see that, in order to renormalize P n (m) to order p, we 
have to solve the first p renormalized equations in the 



hierarchy (jlQ|) . with the second equation solved to order 
p — 1, the third to order p — 2, the p-th to first order. 

Turning to the remnant, P N will contain, order by or- 
der, corrections in the form P R (Z) — P; R (0 an d [Pi(m) — 
Pi(to|Z)] — [Pi (to) — Pi(m|^)] soc , which do not lead to di- 
vergence in Eqs. (|16M17|) and their counterparts for (|10|) . 

Pursuing the analogy with quantum field theory, we see 
that e = plays a role analogous to the upper critical di- 
mension, which suggests us to calculate the renormalized 
probability using an e-expansion approach. We choose 
to keep in [. . .] soc only the scaling part, and, substituting 
Eq. ([II]) into CHI, we get: 



1 n— 1 

AiP„ r (to) = --<z 2 APi(oo) J2 P?(l)(r 



I) 



-l+e 



(18) 



z=i 



For n — > oo, we approximate the sum by an integral; 
defining f(z) = P^(k) with k — n(l — z), we write: 



n-l 

E 

i=i 



P?(l)(m-iy 



dzf(z)[m/n—l+z 



-l+e 



where z min = 1/ra, z max = 1 — 1/n and the factor n e comes 
from n~ 1+e / Az with Az = 1/n the discrete increment in 
the integral. The interesting case is to = n. Integrating 
by parts and expanding in e: 



dzf(z)z- 1+e = e-\n e ~ l)P R (n) + 0(e). 

and then into (fl7jl). we obtain 



Substituting into Eq 
the result: 



and, sending n 



1 



Aq 2 P 1 (oo) 
2e 



K " 1) 



Pi{n) 



P R (n) = Pn" £ + 0(e 2 



(19) 



(20) 



where P = 2e/(Aq 2 ). Notice again the analogy with 
quantum field theory, with P behaving like the fixed 
point value of a renormalized coupling constant. Sub- 
stituting into Eq. ©, we get the stretched exponential 
scaling for the return time PDF: 



P{R q > n) cx exp(-Pn 1 " £ ). 



(21) 



The expectation that the result in [8J] extends to the re- 
turn time statistics is therefore confirmed, and we have 
gained knowledge of the prefactor in the exponent. Con- 
versely, in the range e < 0, no divergence would have 
arisen in Eqs. (|14M15[) . so that no renormalization would 
have been necessary. Hence, the bare perturbation the- 
ory would have been appropriate, with the ground state 
P^ 0) (n) — > Pi(oo) leading to exponential scaling for the 
return time distribution. 

We can use Eqs. (|6I19|) to study the approach to the 
asymptotic regime (|2ip . In this transient regime, it is 
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FIG. 2: Return time probability for e = 0.5 (left) and e = 0.2 
(right, curve a). Curve b on the right refers to exit from the 
initial condition y(0) — ~q, again for e = 0.2. In all cases 
A = 2 and q — 3, corresponding to R q ~ [Pi(oo)] -1 ~ 677. 
Thin lines: numerical integration of y(t) — ~}2 n x n (t), us- 
ing Eq. (J3J) with N = 24. Dotted lines: theory using 
Eqs. (|6I19[) . Heavy lines: theory using Eqs. (|26H27[) h In- 
sert: same results without rescaling, superimposed with the 
stretched exponential fits (circles): 0.004 exp(—0.17?q' 5 ) and 
0.0016 exp(-0.0077?°' 8 ). 



appropriate to use the expression for Pi(n) that is ob- 
tained from Eqs. (|3ll2p . rather than from the asymp- 
totic formula (fT3")) . We see in Fig. [2] that the analytical 
approximation does a pretty good job even up to e = 0.5, 
which is consequence of the fact that the first correction 
to AiP"(n) arises only at O(e). 

Is this stretched exponential scaling? Actually it is 
possible to fit P{R q ) with stretched exponentials, but 
the fit does not match Eq. (|2"Tj) . Using Eq. ([3]), in the 
two cases e = 0.5 and e = 0.2, we evaluate A ~ 0.5 and 
A ~ 0.4; for q = 3, this would correspond to P ~ 0.2 in 
both cases, which goes from a factor of two to two orders 
of magnitude away from the fits in Fig. [2j 

What is happening is that Eq. (fll?)) reaches its asymp- 
totic limit (|2U)) only for n > [P/Pi(oo)] 1 / e , which grows 
very rapidly with q and 1/e. To have an idea, for 
e = 0.2 and q = 3, we would have Pi (00) ~ 0.0015 and 
[P/Pi(oo)] 1 / e ~ 10 10 . Clearly such long return times 
would occur only with vanishingly small probability. 



A. Higher order corrections 

Some idea of the higher orders in the perturbative ex- 
pansion could be obtained studying Eq. (jTTJ) for p = 2. 
We focus on the case m = n and start by analyzing 
whether the second line in Eq. ([T7]). Ai[P ; (n) - Pi(n\l)] 
gives a secular contribution. From Eq. (|16p . we thus 
have to evaluate: 

AlW = Ej=i tf(j)[Pi(n) - Pi(n|j)] soc 
Aifl(n|0=E5=i^010[^(n|0--Pi(n|y)U 



From Eq. (fl~4f and the argument leading to (|15p . 
we see that [Pi{n) — Pi(n\j)] scc ~ \n — and 
[Pi(n\l) - Pi{n\lj)) SDC ~ |i - j|" 1+£ . We know from 
Eq. ® that P, R (j) - j~ e ; we still need P?(j\l). We 
can repeat the steps from Eq. (p~8j) to (| 19[) . substituting 

E^l Pf(j\ l ){ 1 - l) _1+e in the sum of Eq. JTS} and the 
final result is again 

p;(ilo~r e . 

Indicating j = n(l - z) and /(z) = Pf(j),Pf(l\j), the 
leading order behavior for l/n — > of the sums in Eq. 
([2^)1 is therefore: 

Z £ / /(*)[n/J - 1 + z]- 1+t ~ ^- 1+e . 
Jo 

We see that no divergences are present at small I in the 
sum 

n-r 

^PfCOA^W-PKnlO], (23) 
;=i 

so that [Pi(n) — Pi(n\l)] scc = and no renormalizations 
are necessary to this order [see Eq. (fT7|) ]. 

We consider now the remnant, which receives contribu- 
tions at both orders p = 1 and p = 2. The only contribu- 
tion surviving for n — > 00 turns out to be that at order 
p = 2, produced by the sum in Eq. (|23[) . The sum is 
dominated by I — > n; hence: AiP/(n) ~ P^(n) — P\(n) ~ 
— Pi(n), AiP„_i(n|n — 1) ~ — Pi(n|n — 1), and, for suf- 
ficiently large A, Pi(n|n — 1) ^> P\(n). For n — * 00, 
Pi(n|n — 1) = Pi(l) and we estimate: 

P„ N (n) ~ P n R (n)Pr(l). (24) 

Thus, the validity of the renormalized expansion rests 
on the smallness of the exit probability after one step 
A = Pi(l), that behaves like an expansion parameter for 
the theory beside e. In order for the theory to work, it 
is then necessary that the sampling constant A be suffi- 
ciently large. However, this appears to be a rather weak 
constrain, as, already for A = 1, q = 3 and e = 0.5: 
A ~ 0.15. 

Substituting into Eq. (fTT)) for p = 3, we see that 
Pn(n) contributes to A 2 [P;(n) — Pi(n\l)] sec and to renor- 
malization of P„(n), while terms like Pi(lji) contribute 
to P ; N (n|Z) with the same mechanism that lead to Eq. 
(f2"4")l . This suggests that the higher order renormaliza- 
tions to P„(n) are 0(X n ) corrections to prefactors in Eqs. 
(|19H21[) , while the exponent in (f2T))) should remain invari- 
ant. This exponent depends in fact only on the part of 
P;(rc.|k) — P;(n|k)|t ->-oo with the slowest decay, which is 

|io|~ 1+e f°r all p [see discussion leading to Eq. (fTS"]) ]. 

V. PDF STRUCTURE 

The key element of the analysis carried on so far is 
that the conditioned return probability F^(n) ~ Pn~ £ 
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[see Eq. (|20|) ] goes to zero for n — > oo. Recalling the 
discussion at the end of section [II] it is this behavior, in 
contrast to that produced by bare perturbation theory 
P n (n) — > Pi(oo) > 0, that leads to anomalous scaling of 
the return time distribution. This means that the PDF 
Pn(n) determining P^(n) through the relation P^(n) — 
I q / 9 n( ri )dyn, must have vanishing tails at y n > q. The 
fact that for n — > oo p*(n) ^ p(yn), is the signature of 
the long memory of the process. 

An equation for the PDF p n (n) could be derived re- 
peating the steps leading to Eq. Q: 

n-l 

p n {m) = pi(m) + ^2 Pi( l )iPi( m ) ~ MHOL 
i=i 

with obvious definitions for the various PDF's appearing 
in the formula. Taking moments, we obtain equations for 
the conditional mean: 

n—l 

p n (m) = n\{m) + Pi(l)[pi(m) - pi(m\l)}, (25) 
1=1 

and similarly for the higher moments of p n {n). The per- 
turbative analysis of Eq. ([25]) is identical to that of Eq. 
©■ Using Eqs. (|12M13[) . we see that the same pattern of 
divergences is produced: 

[pi(n) - pi(n\l)] scc = -qA\n - 

and this confirms the role of the conditional mean 
p,i(n) — (y n \yo > q) hi the renormalization proce- 
dure. Generalizing Eq. (|25[) to the higher moments 
M n!P (m) = j p n (m)yf rl dy m , in fact, it is possible to see 
that M hp (n) - M hp (n\l) ~ In-Z^ 1 " 6 ) and for e < 0.5, 
the higher moment equations do not have divergent be- 
haviors. 

Using Eq. (f2"5|) , it is possible to renormalize (|2~i|) using 
the same procedure leading from Eq. (TT5)) to (|2"U|) . This 
leads to the result: 

pl(n) = Ml (n) - (q/e)(n< - 1)P„» -> -2/q, (26) 

which is confirmed in Fig. [3] [as in Fig. [5] with Pi(n), 
the expression from Eq. p2[) is adopted here for pi(n)]. 
Notice the constant value for n — > oo of p^(n), much 
larger than the value —qP(y > q) that would be obtained 
subtracting from the equilibrium PDF the values of y 
above threshold. 

Clearly, knowledge of the conditional mean Pn(n) is 
not sufficient by itself to guarantee vanishing PDF tails 
at y > q. However, using such a simple approximation 
for p n (n) as: 

p n (n) ~ (27T)" 1 / 2 exp{-(jM - /i») 2 }, (27) 

with Pn( n ) given by Eq. (|26[) . and substituting into 
Pn{n) = p n (n)dy n and in ([5]), produces results in Fig. 
[2] in some case better than from Eq. ([19]) . One reason for 
this is the better scaling properties of pi(n) — pi(n\l) ~ 
— (Ui(n — /) as compared with Pi(n) — Pi(n\l) (see Fig. 
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FIG. 3: Scaling for the conditional mean p n (n) for q — 3, 
e — 0.2 and A = 2. Dotted line: "bare result" pi(n); heavy 
line: renormalized result p n (n) using Eq. H26[l : thin line: 
numerical simulation using the same parameters of Fig. [2] 



VI. THE NON-GAUSSIAN CASE 

We have seen that the scaling of the return time dis- 
tribution depends on the tails of the conditioned PDF 
p n (n). These in turn are determined by the behavior of a 
bulk quantity like the conditional mean pi(n) = (y n \yo > 
q) . It seems therefore that it is not the extreme statistics 
of the process, rather, its correlation structure that deter- 
mines the return time distribution. It is natural at this 
point to question whether Gaussian statistics is strictly 
necessary for (quasi) stretched exponential scaling. 

We examine the conditions under which pi(oo) — pi(l) 
and therefore also the difference Pi(oo) — Pl(0 in Eq. 
(TTJ}, scale like pi(l) ~ Z~ 1+c . We provide a sufficient 
condition for this scaling in the form of a requirement of 
weak correlation between scales in the process. 

Let us write yo = %o + z an d yi = xi + z, with z the 
result on y of some low-pass filtering at scale I in the 
region in exam. From Eqs. IHU), we have, for large 
I: (z 2 ) ~ l~ 1+e . We can obtain pi(l) and pi(oo) from 
p{y) and p(y, y'; I), which are respectively the equilibrium 
PDF for y and the joint PDF that yi = y and yo = 
y' . Indicating by p < and p > the PDF's for the low-pass 
filtered signal z and for x = y — z, we can write: 

P(y,y',l) = Jdzp < (z)p > (x,x';l\z), , , 
p(y') =fdz P< (z)p > (x'\z). ^ 

We can introduce a function g(y,y', z,l) parameterizing 
the correlation between the small scale components xi = 
yi- z and x = y Q - z: 

p > (x,x';l\z) = [1 + g(y,y',z,l)}p > (x\z)p > (x'\z). 

If correlation between scales are weak and I is large, we 
can Taylor expand p > (x\z) and p{x'\z) around z — 
and consider g(y, y', z, I) a small quantity. With these 



7 



substitutions, Eq. (|28|) becomes: 



p(y,y';i) 
p(v') 



P> (y\0) P> (y'\m- 

P> (y'\0)[l + B(z 2 ) 



(g\y,y';l)+A(z 2 )], 



where A — A{y, y') and B = B(y'). From here, we obtain 
finally: 

p(y\y'; I) ~ p(y)[l + (g\y, y>; I) + (A - B)(z 2 )}. (29) 

The two contributions to pi(l) — Pi(oo) ~ p{y\y'] l)~ p(y) 
are deeply different in nature. The term (A — B)(z 2 ) is 
the direct additive contribution from the long time-scale 
fluctuations, while (g\y, y'\l) probes long time correla- 
tions of the small scale fluctuations. The last contir- 
bution and the direct coupling between fluctuations at 
different scales [l5j|, are typically associated with a mul- 
tifractal structure of the signal 16]. Dominance of the 
additive contribution (z 2 ) ~ l~ 1+e , indicates therefore 
absence of multifractal properties in the signals. The dif- 
ference -Pi(oo) — Pi(l) in Eq- (fl4|) will behave in this 
case as if y were Gaussian, and P(R q ) should become a 
stretched exponential in the large R q limit. 

The fact that fractal objects are not processes with a 
stretched exponential distribution of return times is not a 
surprise. This is easy to see in the case of a middle-third 
Cantor set, in which the return times can be identified 
with the "holes" in the measure: at the n-th generation 
there are 2 n_1 holes of length R q — 3~" and this gives the 
power- law distribution P{R q ) cx R-° with D = In 2/ In 3 
the fractal dimension of the set. Actually, there have 
been some recent attempts to characterize multifractal 
sets by a return time spectrum, beside the more standard 
singularity and dimension spectra [l7l [l8| . 

The prediction that the return time statistics is domi- 
nated by the correlation structure of the process is con- 
firmed by numerical simulation of non-Gaussian power- 
law correlated processes, as illustrated in Fig. 0J The 
four curves in figure are all characterized by the same 
power-law correlation with e = 0.2, but are generated 
with different forms of intermittency in the sub-processes 
x n (t) in y(t) — J2 n x n(t). The intermittency is generated 
letting the noise amplitude in the Langevin equation gov- 
erning each process fluctuate on the time scale of the 
process: 

X n (t) = -TnXnit)- 1 + b n (t)£(t), (mm) = S(t) 

with (6^) ~ 2a 2 ri /T n . In the multifractal case, the noise 
amplitude fluctuations are generated by a multiplicative 
process of the kind utilized in [16j. In the other cases, 
the noise amplitude fluctuations are independent. In the 
first case, the fluctuations are tuned to produce the same 
kurtosis (y A ) ~ 12.5 as in the fractal case. In the second 
case, the intermittency grows with scale [which would 
produce non-trivial scaling of the higher diffusion expo- 
nents for r(t) = J y(r)dT]. In the third case, the signal 
is Gaussian. As expected, the return time distributions 
of the non-multifractal signals collapse on one another if 
one rescales with R q , while the multifractal one leads to 
a distribution that is closer to a power law. 
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FIG. 4: Return time probability for different non-Gaussian 
processes. In all cases, e = 0.2 and A = 1. Case a: three non- 
multifractal signals with different types of non-Gaussianity; 
case b: multifractal signal. 



VII. PERMANENCE TIME DISTRIBUTION 

Let us conclude the analysis turning to the permanence 
times and verifying that their distribution is character- 
ized by stretched exponential scaling, as extension of the 
theory in Q would suggest. The analysis is limited to 
the Gaussian case. The analog of Eq. ([5]) in the case of 
the permanence times S q reads: 

P(S q >t) = P-\ V > q)P(y(T) >q,re [0,i]), 

where the factor P~ 1 (y > q) gives the condition that the 
stochastic variable is initially above threshold. Let us 
isolate in y(r) its average in [0, t]: y(r) = z + x(t), where 
z = t^ 1 J Q * y(r)dr. We then obtain: 

P(S q >t) = p- 1 (y>q) 
x dz P< (z)P(x(T) >q-z,re [0,t]\z). W 

The PDF for z, for large t, is obtained eliminating fre- 
quencies \u>\ > t~ x in the power spectrum for y; from Eq. 
©: (z 2 ) - f, t|<1 C w dw - t~ 1+e . In the same limit, 

the condition on z in P(x(t) > q — z,t G [0, t]|z) can be 
disregarded for Gaussian statistics. 

To evaluate Eq. (|30|) , we consider the simpler problem 
of discrete sampling in time: t — > t n = nA. This allows 
us to write: 

1 > P{x(t) > q - z, t G [0, t]) > [P{x >q- z)] t/A . 

Substituting into Eq. (f30|) , the first inequality allows us 
to write, to leading order in q and t: 



P(S q > t) < exp(-KqH L - e ), 



(31) 



where we have estimated p<(z) ~ exp(— Kz 2 t x e ) and 
then, for large q: P(z > q) ~ exp(— Kq 2 t 1 ~ e ). 
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FIG. 5: Permanence time distribution from numerical sim- 
ulation of a power-law correlated time series (thin lines) and 
stretched exponential fit (heavy lines) for q = 3 and three dif- 
ferent values of e: (a) e — 0.8; (b) e = 0.4; (c) e = 0.2. Insert: 
scaling of {y M Ax\S q ) vs. S q for e = 0.4. 



Passing to the second inequality, making the substitu- 
tion in Eq. ijSDJ): P(x(t) > q- z,t G [0,t]) -> [P(x > 
q — z)]*/ A , the integrand in that formula will take the 
form: 

exp{-Kz 2 t 1 - £ + (t/A) ln[l -P(x <q- z)]}. 

For t/A — > oo, the contribution to the integral will be 
from values of z for which ln[l — P(x < q — z) is small, 
i.e. z — q is large; we can approximate: 

ln[l - P(x < q - zj] exp(-K(q - z) 2 t 1_e ). 

The integrand in Eq. (|30|) will then take the form 

exp{-Kz 2 t 1 - £ - K'texp(-K(q - z) 2 i 1_£ )}, 

which is peaked, for t — > oo, at z ~ q+\fe\nt. Estimating 
the integral in Eq. (f30|) by steepest descent gives the 
result: exp(— Ket 1 ~ t lni). Combining with Eq. (|3"Tj) . we 
obtain the bound, valid to leading order in q and t: 

cxpi-Ket^niit) < P{S q > t) < expf-Jf^ 2 * 1-6 ), 

(32) 

which is similar in form to the one in [8j]. Contrary to 
the case of the return times, the value of the sampling 
constant A is not crucial to the theory. This is confirmed 
by the fact that, in the present limit, all dependence on 
A, accounted for by K', disappears. We compare in Fig. 
[5] with the result of numerical simulation using Eq. (|3|) 
and see that stretched exponential scaling is compatible 
with the permanence time distributions in the range con- 
sidered. 

Notice that replacing the upper bound in Eq. (|32|) with 
equality would imply: 

P(S q > t) ~ P(Z >q)~ P(y > qt 1 ^), 



in other words, the probability of a single peak of height 
qSq~ e would be the same as that of a permanence S q . 
This is not surprising: given an initial condition y(0) 3> 
g> 1, from Eqs. (|11I13[) . p(y(t)\y(Q)) would be narrowly 
peaked around y(0)At~ 1+e and the time it takes to y(t) 
to go below q would be t ~ (^4y(0)/g)~. Substituting 
y(0) = qSq~ e , we obtain precisely t ~ S q . As confirmed 
also in Fig. one expects therefore that longer perma- 
nences above threshold be associated with higher peaks. 

VIII. CONCLUSION 

The analysis carried on this paper confirms the obser- 
vation in 0, [5], Q (among others) that long correlations 
in stochastic processes lead to return time distributions 
with scaling close to stretched exponential. Similar prop- 
erties, albeit with different mechanisms, are confirmed for 
the permanence time distribution [see Eq. (|32p . 

These results are consistent with the extension to 
thresholds different from zero, of the bounds derived in 
Q on no zero-crossing probabilities. The same analysis 
also suggests that what can be observed in experimen- 
tal time series is only a very slow transient regime. The 
stretched exponential scaling predicted in [8j is achieved 
only as an asymptotic limit, requiring exceedingly long 
return times and vanishing probabilities. However, as 
the theory is based on an e-expansion, it is not ruled out 
that the asymptotic limit may occur earlier in the range 
e — * 1, where new physics may become important. 

From the practical point of view, it is probably irrele- 
vant whether a return time distribution that can be fitted 
by a stretched exponential is really a stretched exponen- 
tial. More important, the present theory provides ap- 
proximate expressions for the return time distribution, 
valid for any value of the argument and working well up 
to e ~ 0.5, i.e. the middle of the range considered [see 
Eqs. © and or (|26I27|) ]. The theory is limited to 
discrete sampling and is actually in the form of a double 
expansion in e and the parameter A = P(yi > q\yo > q), 
which is heavily dependent on the sampling constant A 
[see Eq. (f24|)]. However, for large enough q, in the range 
of e in which the theory works, considering A in the scal- 
ing range for the correlation C(t) appears to be sufficient. 

Another important fact is that, although the approx- 
imate expressions in the present paper are derived in 
the Gaussian case, they continue to be valid for non- 
Gaussian processes, provided one rescales the relevant 
quantities by the mean return time R q : R q — > R q /R q 
and P{R q ) — * P(R q )R q . Basically, only multifractal pro- 
cesses are excluded. (The importance in this context of 
rescaling by R q was first pointed out in 3). 

The mechanism for the non-exponential scaling of the 
return time distribution appears to be that trajectories 
originating from an above-threshold event are distributed 
around a mean (y(t)\y(0) > q) that decays very slowly 
with time. When imposing that the trajectories remain 
below threshold up to the return time R qi loosely speak- 
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ing, this slow decay produces correlations between the 
conditions of no-exit at different times that cannot be 
treated as independent. The important point is that, as 
long as no slower scaling quantities are characterizing the 
process (due e.g. to multifractality) , this mechanism will 
be insensitive to whether or not the process is Gaussian 
[see Eq. (j2"9")l and following discussion]. This insensi- 
tivity on the tail structure of the statistics was recently 
observed in fl9j | . 

From the conceptual point of view, it is interesting 
that the cumulative effect of the correlation between the 
below-threshold conditions becomes manifest through 
secular behaviors in the equations for the evolution of 
the trajectory distributions. It is also interesting that 
the most natural way to treat these behaviors is renor- 
malization, with a strategy similar to [141 ]. In particular, 
the stretched exponential scaling limit appears to be as- 



sociated with the fixed point of the renormalized theory. 
The range e < 0, corresponding to a process with finite 
correlation time and exponential scaling for the return 
distribution, conversely, is associated with a trivial the- 
ory that does not require renormalization. 

It is to be noticed that equations similar to the ones 
considered in the present paper arise in the context of 
avalanche models [201 ]. where they have been treated as 
well within an e-expansion approach. 
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